A magnetohydrodynamics theory was used to study the dynamic response of the neutral atmosphere to a geomagnetic storm. A full set of magnetohydrodynamic equations appropriate for the present problem is derived and their various orders of approximation are discussed in some detail. In order to demonstrate the usefulness of this theoretical model, the May 1967 geomagnetic storm data were used in the resulting set of nonlinear, time dependent, partial differential magnetohydrodynamic equations to calculate variations of the thermosphere due to the storm. The numerical results are presented for wind speeds, electric field strength, and amount of joule heating at a constant altitude for the data recorded. Data show that the strongest thermospheric responses are at the polar region becoming weaker in the equatorial region. This may lead to the speculation that a thermospheric wave is generated in the polar region due to the geomagnetic storm which propagates towards the equator.